function [fea,sco] = mrmr_miq_d(d, f, K)
% function [fea] = mrmr_miq_d(d, f, K)
%
% MIQ scheme according to MRMR
% Mutual Information Quotient
%
% By Hanchuan Peng
% April 16, 2003
% 
% TODO:
% Output `sco` (score) is of problem right now. Need revision later.
% Note that there is `KMAX` (max allowed size of feature sub-set) in the
% code.
% Revised by raywayzhao, 2011

% bdisp=0;

nd = size(d,2);
% nc = size(d,1);

% t1=cputime;
for i=1:nd
    t(i) = mutualinfo(d(:,i), f);
end
% fprintf('calculate the marginal dmi costs %5.1fs.\n', cputime-t1);

[tmp_score, idxs] = sort(-t);
fea_base = idxs(1:K);
fea(1) = idxs(1);
KMAX = min(1000,nd); %500 %20000
sco = nan(1,KMAX);
sco(idxs(1)) = -tmp_score(1);
idxleft = idxs(2:KMAX);

% k=1;
% if bdisp==1
%     fprintf('k=1 cost_time=(N/A) cur_fea=%d #left_cand=%d\n', ...
%         fea(k), length(idxleft));
% end

for k=2:K
%     t1=cputime;
    ncand = length(idxleft);
    curlastfea = length(fea);
    
    %TODO: Consider pre-allocating `t_mi`, `mi_array`, `c_mi` here
    for i=1:ncand
        t_mi(i) = mutualinfo(d(:,idxleft(i)), f);
        mi_array(idxleft(i),curlastfea) = getmultimi(d(:,fea(curlastfea)), d(:,idxleft(i)));
        c_mi(i) = mean(mi_array(idxleft(i), :));
    end
    
    %   [tmp_score, fea(k)] = max(t_mi(1:ncand) ./ c_mi(1:ncand));
    [tmp_score, fea(k)] = max(t_mi(1:ncand) ./ (c_mi(1:ncand) + 0.01));
    tmpidx = fea(k); fea(k) = idxleft(tmpidx); idxleft(tmpidx) = [];
    sco(fea(k)) = tmp_score(1);
    
%     if bdisp==1
%         fprintf('k=%d cost_time=%5.4f cur_fea=%d #left_cand=%d\n', ...
%             k, cputime-t1, fea(k), length(idxleft));
%     end
end

end

%=====================================
function c = getmultimi(da, dt)

for i=1:size(da,2)
    c(i) = mutualinfo(da(:,i), dt);
end

end
